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Predictions for Triple Stars with and without a Pulsar in 
Star Clusters 
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ABSTRACT 

Though about 80 pulsar binaries have been detected in globular clusters so far, no 
pulsar has been found in a triple system in which all three objects are of comparable 
mass. Here we present predictions for the abundance of such triple systems, and for the 
most likely characteristics of these systems. Our predictions are based on an extensive 
set of more than 500 direct simulations of star clusters with primordial binaries, and 
a number of additional runs containing primordial triples. Our simulations employ a 
number N to t of equal mass stars from N to t = 512 to N to t = 19661 and a primordial 
binary fraction from — 50%. In addition, we validate our results against simulations 
with N — 19661 that include a mass spectrum with a turn-off mass at 0.8 Mq, 
appropriate to describe the old stellar populations of galactic globular clusters. Based 
on our simulations, we expect that typical triple abundances in the core of a dense 
cluster are two orders of magnitude lower than the binary abundances, which in itself 
already suggests that we don't have to wait too long for the first comparable-mass 
triple with a pulsar to be detected. 

Key words: globular clusters - pulsars, general - methods: N-body simulations - 
stellar dynamics. 



1 INTRODUCTION 

Numerical simulations of star clusters routinely pro- 
duce dynamically formed triple s tars, especially i n the 
presence of primordial binaries llMcMillan et al.l Il99d ; 
iHeggie fc Aarsethlll992l ; lMcMillan fc Hutlll994 ). Given that 
most globular clusters have a significant population of pri- 
mordial binaries, we would expect some pulsar binaries in 
globular clusters to be actually part of a triple system. In 
such a case, most likely all stars have comparable masses, 
given the fact that encounters between binaries, and also be- 
tween binaries and single stars, are likely to expel the lightest 
star, and hence subsequent encounters tend to concentrate 
the most massive stars in binaries and triples. Such a system 
we refer to as a pulsar triple, i.e. a triple system in which 
all components have comparable mass and one is a pulsar. 
(This phrase is preferable to the alternative triple pulsar, 
which would strongly suggest a triple system in which all 
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three components are pulsars.). The pulsar that is detected 
through its radio emission can then either be one of the in- 
ner two stars, or it can be the third star in orbit around the 
other two. The possible set of configurations is quite large: 
each of the other two stars can be a neutron star, a massive 
white dwarf, possibly a main sequence star, and perhaps 
even a black hole. So far, no pulsar triple with compara- 
ble masses has been discovered, among the more than 120 
pulsars found in globular clusters. 

Ideally one would like to explore the formation of pulsar 
triples using detailed N-body simulations which include di- 
rect integration of the orbits and stellar evolution (therefore 
using one particle per star). Unfortunately this approach is 
not computationally feasible today nor it will be in the near 
future. In fact, even if we consider a pulsar rich globular 
cluster such a s Terza n 5, the population synthesis models of 
ivanova et all il2007ft , which include a simplified treatment 
of the dynamics of the system, predict a total of « 120 
neutron-star binaries (that is binaries with at least one neu- 
tron star) and thus about one neutron-star triple, consider- 
ing that stable triples are about 1% of stable binaries (e.g. 
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see lAarsethll200ll ; see also Sec. [6}. This immediately high- 
lights that, in order to gather enough statistics, one would 
need to run several tens of simulations of star clusters such 
as Terzan 5, which has a present mass of 3.7- 10 5 M Q and a 
number of stars N > 4 • 10 s . This is outside the capabilities 
of current hardware, including the GRAPE6: the largest di- 
rect N-body simulations c arried out to date for at least a re- 
laxation time are tho se by Baumgardt fc Ma kino ( 2003]) and 
have N=131072 (see iHut fe TrentfeOOTT l. Thus one is forced 
to introduce some simplications to address the formation of 
triples stars with a pulsar. As a first step we use equal- 
masses calculations as discussed in Sec. |3] These provide 
a reasonable approximation for the dynamics of old globu- 
lar clusters, where the turn-off mass is below IMq and the 
typical mass ratio between two stars in the cluster is around 
2:1. Such an approach appears to be preferable to the use 
of simulations that start with very massive stars and up to 
several thousands particles, as is instea d appropriate for the 
study of young open clu sters (e.g. see Ide la Fuente Marcos! 
ll99rj ; lHurlevet al.ll2005l) . In addition, we validate our equal- 
masses results against simulations with a mass spectrum 
that has a turn-off mass at 0.8Mq, consistent with the age 
of stellar populations in galactic globular clusters. 

The aim of this paper is two-fold. First, we characterize 
the typical fraction of triples that are formed dynamically 
in a star cluster with a sizeable population of primordial 
binaries, comparing their number density with that mea- 
sured at late times in simulations that start with primordial 
triples. Second, we take advantage of the fact that triple sys- 
tems with one or more pulsars have most likely a dynamical 
and not primordial origin. We therefore apply our results to 
make an estimate of the likely time we have to wait until 
the first pulsar triple will be found as well as to predict its 
most likely characteristics. 

The paper is organized as follows. In the next two 
sections we present an overview of observations and sim- 
ulations, respectively. In Sec. 4 we present the initial 
conditions for the simulations presented in this paper. 
Sec. 5 briefly summarizes the global evolution of a star 
cluster with primord i al bin a ries, extensively discussed i n 



measured to 10— 100/is), if they were part of a triple sys- 
tem, the other two components would be easily detected, 
even if one o f the components w as low-mass or in a long- 
period orbit dThorsett et al.lll999l '). 

Since binary MSPs typically have orbital periods of sev- 
eral hours to several days, a pulsar triple system with the 
MSP in the internal orbit would likely be initially identi- 
fied as a normal binary pulsar system. Timing observations 
would expose the third body in the external orbit after a 
time equal to a few percent of the external orbital period, 
which is likely to be less than a year. 

Alternatively, a pulsar triple system could have the 
MSP in a relatively long-period (months to years) external 
orbit around a compact inner binary likely comprised of neu- 
tron stars, massive white dwarfs, or a combination of both, 
and possibly main-sequence stars. The MSP in such a sys- 
tem would be initially identified as a long-period binary or 
possibly even an isolated pulsar. Timing observations would 
relatively quickly reveal the MSP's external orbit, but un- 
less strong orbital perturbations due to the inner binary are 
present, the internal orbit might never be detected and the 
MSP "companion" will be assumed to be a single massive 
star or compact object instead of a binary. 

It is possible that stellar interactions could produce a 
triple system where two of the members are radio MSPs. The 
MSP in the external orbit would be relatively easy to detect 
due to the likely small orbital accelerations present over the 
short time intervals (i.e. hours) used for pulsar searches. For 
an MSP in the compact inner orbit, though, strong orbital 
accelerations would make the initial detection of the MSP 
very difficult, likely requiring specializ ed algorithms and e x- 
tremely large amounts of computing l|Ransom et al.ll2003h . 

Finally one could ask what is the a-priori chance to find 
a triple system with two of the three bodies being pulsar. 
We expect that this probability is significantly smaller than 
that of finding a triple system with just one pulsar. In fact 
among the ~100 binaries containing a pulsar, none is a dou- 
ble pulsar. This suggests that the probability for a triple 
with two pulsars is at least two orders of magnitude smaller 
than the probability of finding a single pulsar triple. The 
Heggie. Trenti fc HutJ (|2006l ); iTrenti. Heggie fc Hutj l|2007h ; expected number of triples with two pulsars {Nt_ pp ) can be 



Trenti et al l (|2007h . In Sec. 6 we develop a simple model estimated as: 



to describe the formation and destruction of triples. We 
use the model to interpret the observed triples fraction in 
our runs, attempting to characterize its A^-dependence. In 
Sec. 7 we discuss in detail the properties of these dynam- 
ically formed triples. Sec. 8 describes the results of some 
additional runs where we have started with the presence of 
primordial triples. Sec. 9 presents the prospects for detection 
of a triple system with a pulsar and Sec 10 sums up. 



2 OBSERVATIONS 

There are currently ~130 pulsars known in globular clus- 
ters, the majori ty of which are in bin ary systems (for a re- 
cent review, see ICamilo fc Rasiol E^OOsI 1 ! Since almost all of 
these systems are millisecond pulsars (MSPs) with extraor- 
dinary timing precisions (pulse arrival times are typically 



N p N B . 



1 See P. Freire's cluster pulsar catalog 
http : //www2 .naic . edu/~pf reire/GCpsr .html 



at 



where N p is the number of pulsars, N 3 is the number of single 
stars, Nb is the number of binaries, Nb. P s is the number of 
binaries with one pulsar, Nt is the total number of triples of 
the system and 77 is the fraction of pulsars not bound to an y 
companion in globular clusters (77 ~ 0.1 from lFreirell2005h . 
If we assume for a globular cluster Np m 10, N B « 10 5 , 
N B . pa » 10, N B ~ 10 4 , Nt ~ 10 2 we obtain N T . PP » 10 -6 . 
Therefore it is unlikely that any double-pulsar triple will be 
detected in globular clusters. 



3 SIMULATIONS 

The dense cores of globular clusters represent natural lab- 
oratories for studying exotic stellar populations, frequently 
living in multiple systems, such as blue stragglers, low mass 
X-ray binaries, cataclysmic variables, millisecond pulsars, 
stellar and, possibly, intermediate mass black holes. Most of 
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these objects form through the combined effects of stellar 
evolution and stellar dynamics, especially due to interaction 
and evolution of binary stars, that may constit ute up to 50% 
of th e core mass of a globula r cluster (e.g., see Albrow et al.l 
l200ll ; iBellazzini et al] |2002| ; iPulone et atl |2003| ). Observa- 
tional evi dence for triple stars is also s t arting to accumulate 
(e.g., see iRaboud fe Mermilliodl Il998l ; ISterzik fc Tokovininl 
120021 ). 

From the theoretical point of view, a detailed model- 
ing of dense stellar systems is extremely challenging due to 
the huge dynamical range involved. Even neglecting stellar 
evolution and hydrodynamics and focusing on gravitational 
interactions only, the local orbital timescale of a binary star 
is several orders of magnitude smaller than the global re- 
laxation timescale (hard binaries have an orbital period of 
a few hours, while the half mass relaxation time may be up 
to a few billion years) . 

Numerical simulations with primordial binaries have 
thus been performed either using approximate algorithms 
such as Fokker Planck or Mont e Carlo methods (iGao et all 
Il99ll : iGiersz fc Spurzeml l2000l : iFreeeau et aD I2003T I. that 
have to rely o n physical inputs on th e in teraction cross sec- 
tions (but see lFregeau fc Rasioll2007l and lGiersz fc Spurzeml 
120031 for Fokker Planck codes that treats binary interactions 
via direct N-Body integration) , or using direct N-Body codes 
but employin g, until recently, only a limited number of parti- 
cles N as 10 3 (McMillan et al.iri990l: iHeggie fc Aarseth1(l992l : 



iMcMillan fc Hutlll994j : Ide la Fuente Marcos! 1 1996h . Studies 
of hierarchical systems in open cluste rs have been carried out 
bv lde la Fuente Marcos et alJ (|l997h . A few direct numerical 
simulations with N fts 10 4 aimed at studying the formation 
and evolution evoluti on of hierarchical sys t ems have been 
recent ly perfo rmed by Aarseth fc Mardlma (|2000h ; lAarsethl 
|200ll ) and by lAarsethl (|2000l . l2004h ~ 

To enhance the statistical base and remove some of 
the limitations of previous d irect N-body investigations , 
we have started a projec t (|Heggie. Trenti fc Hut 120061 : 
iTrenti. Heggie fc HutJ 120071. herea fter HTH,THH respec- 
tively; see also ITrenti et all l2007h to study the evolution 
of star clusters with primordial binaries by means of sev- 
eral hundred direct numerical simulations employing up to 
Ntot = 19661 particles. Our aim is to provide a milestone for 
comparison and validations of approximate methods and to 
investigate in detail the complex dynamical interactions of 
multiple stellar systems, that cannot so easily be character- 
ized in terms of approximate treatments like Fokker Planck 
or Monte Carlo methods. 

In this paper we take advantage of the extensive set of 
numerical simulations (more than 500) that we have per- 
formed (see HTH,THH) and we estimate the frequency and 
the orbital properties of dynamically formed triple stars. 
We discuss their dependence on the dynamical state of the 
star cluster (pre and post-collapse), on the initial binary 
ratio and on the number of particles TV. We complement 
these simulations (i) with previously unpublished runs that 
include a mass spectrum appropriate to represent the colli- 
sional evolution of old stellar populations in globular clus- 
ters (that is with a turn-off mass below IMq) and (ii) with a 
few equal mass runs with a significant population of primor- 
dial triples. Although most of our simulations are limited 
to equal mass stars without stellar evolution, we briefly dis- 
cuss in the concluding section the applicability of our results 



to the prediction of the properties of globular cluster triple 
systems of comparable mass that include one, or more, mil- 
lisecond pulsars. 



4 INITIAL CONDITIONS 

The general properties of the numerical simulations consid- 
ered in this paper to study the properties of triple systems 
have been presented in detail in HTH a nd THH. To su mma- 
rize, we use Aarseth's NBODY-6 code (lAarsethll2003l ) con- 
sidering stars of equal mass and no stellar evolution. The 
initial distribution is either a Plummer model (HTH; in this 
case the system is considered isolated) or King models with 
concentration parameter Wo = 3, 7, 11 (THH; here the ef- 
fects of a self-consistent galactic tidal field are also taken into 
account). The number of objects used ranges from 256 to 
16384; here N = N s + N b , where N b and N s are the initial 
number of binaries and single stars, respectively. The total 
number of stars is Ntot — N s + 2N b . In our standard runs 
(see HTH and THH), the primordial binary population has 
a fractional abundance (/ = N b /(N S + N b )) from to 50%, 
with an interal binding energy distribution flat in log scale 
in the range sa 5 fcT-680 kT, where (3/2)feT is the mean 
kinetic energy per particle of the system (the binaries being 
replaced by their barycenter). In addition we have performed 
a few simulations with N = 8192 and / = 20% that start 
from a King Wo = 7 model (with a self consistent tidal field) 
and adopt an extended binding energy range (from 5 kT to 
10 4 kT) for primordial binaries. The aim of these runs is 
to specifically investigate the formation of triples originated 
from neutron star binaries with a degenerate companion. 

In addition to these equal mass runs we also consider 
two tidally limited runs (Wo = 7 and Wo = 11 King models) 
with a mass spectrum. For these simulations, which have 
N = 16384 and / = 10%, the individual particle mass is 
drawn from an initial mass function appropriate to study 
the late evolutionary stages of star clusters, when stars are 
~ 10 Gyr old. To bu ild the initial mass fu nction for the 
run we first consider a lMiller fc Scald (|l979h mass function 
in the interval [0.2, 10]A/q, and then we proceed to take 
into account the effect of stellar evolution for massive stars 
by reducing the mass of particles above the turn-off (as- 

-ff = 0.8 . 
of lHurlev eta!] <|2000h . 

All our results are presented using standard units 
|Heggie fc Mathieu|[l98rj) in which 



su med at mtoff = 0.8 Mq) accordingly to the prescription 



G = M 



-4E T 



where G is the gravitational constant, M the total mass, 
and Et the total energy of the system of bound objects. In 
other words, Et does not include the internal binding energy 
of the binaries, only the kinetic energy of their center-of- 
mass motion and the potential energy contribution where 
each binary is considered to be a point mass. The cor- 
responding unit of time is t d = GM 5/2 /(-4E T ) 3/2 = 1 
l|Heggie fc Mathieu|[l98r3 ). For the relaxation time, we use 
the following expression (|Spitzerlll987h 



0.138/Vr 



3/2 



ln(0.11/V) ' 

We recall that as our simulations consider gravitational in- 
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teractions only, our results, expressed in terms of the dimen- 
sionless units described above, can be applied to any physi- 
cal choice for the total mass and scale radius of the system. 
The dependence of the triples' properties on the number of 
particles employed is discussed in Sec.|6]and provides a way 
to extrapolate the results of our simulations to a number of 
particles realistic for standard globular clusters. 

For reference in physical units, a globular cluster (de- 
scribed by a Plummer model) with N = 3 • 10 stars, a total 
mass of M = 3 • 10 5 Mq and a half- mass radius of 4 pc has a 
half-mass relaxation time t r h ~ 8.5 • 10 8 yr\ in this cluster a 
binary, formed by equal mass stars each of mass 1 Mq , with 
binding energy of 1 kT has a semi-major axis of w 10 AU 
and an orbital period po ~ 20 yrs. 



5 GLOBAL EVOLUTION 

The evolution of the system in our runs is driven by the 
balance of two competing phenomena: the tendency of the 
core to contract (undergoing a gravothermal collapse) and 
the generation of energy due to binary-binary and binary- 
single interactions, that eventually halts the core contraction 
and fuels the half mass radius expansion. We can broadly 
identify two phases, extensively described in HTH and THH: 
an initial core radius adjustment and a quasi-steady binary 
burning phase. 

The first "adjustment" transient is characterized by the 
evolution of the core toward a quasi-equilibrium configura- 
tion that will be maintained during the subsequent binary 
burning phase. In this phase we may have, depending on 
the initial conditions, a contraction (e.g., starting from a 
Plummer model, or from King models with Wo < 7) or an 
expan sion, if the core is initially too small (|Fregeau et al.l 
120031 and the discussion of the Wo = H runs in THH). 
For equal mass stars this transient lasts up to ~ 10 t r h(0) 
and the expected core radius value at the end of this initial 
adjustment is well modeled in terms of ge neral theoretical 
considerations (Ves perini fc Chernofd[l994r ). Essentially the 
efficiency of the production of energy due to binary burning 
depends on the core density: if the core density is too low, 
the energy production is inefficient and the core shrinks; on 
the other hand a very high core density leads to an excessive 
energy generation, with a consequent core expansion. 

After the "adjustment" transient, the details of the ini- 
tial conditions are largely erased from the system and a self 
similar evolution sets in: the cluster expands by keeping the 
ratio of the core to half mass radius almost constant. The 
fuel for the expansion is provided by the hardening and de- 
struction of the primordial binary population, that can sus- 
tain this phase for about 100 t r h(0) in isolated models (thus 
for a time much longer than the age of the universe for a typ- 
ical globular cluster) if the initial binary fraction is above 
10%. 

As a result of four body interactions, basically happen- 
ing in the core of the system only, where the interaction 
probability is highest, a (small) fraction of stable triple stars 
is formed. In the next Section we discuss the properties of 
these multiple systems. 



6 TRIPLE FORMATION AND ABUNDANCES 

Our runs start without primordial triples, but in a few 
relaxation times stable triples (identified as stable in 
NBODY6 using a semi-ana lytic approach based on th e bi- 
nary tides problem - see lAarseth fc Mardlind |2000| and 
iMardling fc AarsethllioOlh are formed via dynamical inter- 
actions in binary-binary encounters following ejection of one 
star. We recall that stable triples cannot be formed due 
to th ree body encounters (see, for example, iHeggie fc Hutl 
l2003h . 

After a first rapid rise of the number of triples that 
happens in the first few relaxation times, a quasi stationary 
regime sets in and the number of triples is approximately 
proportional to the number of binaries. This regime lasts 
for about 30 relaxation times, a time longer than the age 
of the universe for a typical globular cluster. Eventually, in 
isolated runs, as the evolution proceeds and the reservoir of 
primordial binaries is being depleted, the number of triples 
also drops. Fig [T] illustrates the evolution of the number of 
triples in an isolated run (N=16384, 10% primordial bina- 
ries, starting from a Plummer model) up to w 90 t r h(0)- The 
drop in N t /N for t > 30t r h(0) is clear, while the fraction of 
triples to binaries declines more slowly. Fig. [2] shows instead 
the number of triples for a run that includes a galactic tidal 
field, starting from a King Wo = 7 model with 20% pri- 
mordial binaries and N= 16384. The tidal field destroys the 
cluster in about 30 t r h{0), thus even in the last stages of the 
evolution the binary fraction remains high (as single stars 
preferentially escape from the cluster) and so does the triple 
fraction. The asymptotic fraction of triple to single stars for 
our runs with equal-mass particles is given in Tab. \T\ This 
has been obtained by averaging the triple to single stars ra- 
tio after the core contraption (t > t cc ) and, in the case of 
tidal field runs, stopping when the number of stars in the 
system is reduced below 10% of its initial value. 

The evolution of the system is very similar even when 
a mass spectrum is introduced in the simulations (see Fig. 
and Tab. Within the statistical uncertainties associated 
to the small number of triples present at any time in the 
system (N t < 10), the behavior of these runs is completely 
consistent with that observed in their equal-mass particles 
counterparts. These experiments with a mass spectrum con- 
firm that equal mass runs are a very reasonable approxima- 
tionas a first step for characterizing the properties of pulsar 
triples in systems where the turn-off mass is below IMq. 

The evolution of the number of triples can be under- 
stood in terms of the balance between the formation (r/) 
and the destruction (rd) rate for triples. In first approxima- 
tion we can write locally: 

r f = a f -p 2 b , (2) 

rd = Pt(ctdB ■ pb + ctds ■ Ps), (3) 

where p s , pb, pt are the density of singles, binaries and triples 
respectively, while the m are numerical coefficients that de- 
pends on the cross section for the process considered. 

With the aid of Eqs. (J2J3J) we can explain the behavior 
of Fig. [2] During the initial core contraction phase the for- 
mation rate 77 is increasing as binaries accumulate in the 
core, while the destruction rate is smaller due to the rel- 
ative absence of triples, whose numbers therefore increase. 
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Figure 1. Evolution of the number of triples with respect to the 
number of particles TV (upper panel), to the number of binaries 
(middle panel) and number of binaries to the number of particles 
(lower panel) for a simulation with TV = 16384 starting from 
a Plummer model (isolated, no tidal field) and 10% primordial 
binaries. Due to the absence of a tidal field the mass loss timescale 
is significantly longer than in the run presented in Fig. [2] and the 
drop in the number of triples at later stages of the evolution is 
apparent. See HTH for more details on the run. 
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Figure 2. Evolution of the number of triples like in Fig. [T]but for 
a simulation with TV = 16384 (with a tidal field) starting from a 
King Wo = 7 model and 20% primordial binaries. The time is in 
units of the initial half mass relaxation time t r f,(0). After about 
30t r h(0) the cluster is completely dissolved by the tidal field. See 
THH for more details on the run. 



Figure 3. Evolution of the number of triples like in Fig.[T]but for 
two simulations with TV = 16384 (with a tidal field) starting from 
a King Wo = 7 model (red solid line) and Wo = 11 model (blue 
dashed line) with 10% primo rdial binaries. T hese simulations have 
particle masses drawn from a lMiller"fc Scald jl97Sh mass function 
evolved to have a turn-off mass at 0.8Mq. 



Around core collapse the formation rate reaches its peak and 
in the mean time an equilibrium regime sets in where the 
density of triples can be obtained by equating the formation 
and destruction rate (77 = r<j): 



Pi 



a f 



OldB ■ Pb + OLdS ■ Ps 



(4) 



We can further assume, based on geometrical arguments, 
that the cross section for a binary-triple encounter is greater 
than that of a triple-single encounter. In addition we recall 
that after core collapse in simulations starting with / < 10% 
the central binary density at least equals that of single stars 
(e.g. see Fig. 21 in HTH). This leads to a simplification of 
Eq. ©: 



Pi 



a f 



Pi 



OldB ■ Pb 



OldB 



Pb. 



(•») 



This theoretical prediction is nicely confirmed in Fig. U 
where the abundance of triples after core collapse for runs 
starting from TV = 4096 is shown as a function of the frac- 
tion of primordial binaries: here the structural properties of 
the system are nearly identical for / > 10% (see Figs. 17 & 
18 in HTH) and there is a linear dependence of TV t on /. 

Finally at later stages of the evolution of the system 
(i.e. well after the end of core-contraction phase), pb in the 
core drops as binaries are being depleted, and the number 
of triples also decreases. 

The fractional abundance of triples as a function of the 
number of particles used, for runs starting with / = 10%, 
either from a Plummer or from a King Wo = 7 model, is 
depicted in Figs. 15181 Especially for low N runs, before the 
core collapse the triple abundance is lower than after it. 
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Figure 4. Fractional abundance of triples averaged after the core 
contraction for isolated runs starting from a Plummer model and 
different primordial binary ratios. Each dot represents a different 
run, while the squares, with a la error bar, are the average values 
of the fractional abundances of triples associated with a given 
primordial binary ratio /. The number of triples appears to be 
linearly proportional to the primordial binary fraction. 



This is because the time for core collapse (in units of the 
relaxation time) increases with N (see HTH, Fig. 17), so 
that when N is low enough the system has not yet reached 
a balance between the formation and destruction rate of 
triples. In the runs starting from a King model we observe 
a marginally higher ratio of triples. This is due to the fact 
that the runs starting from a King model have a tidal field 
(see THH), so that the system steadily loses stars in its 
outskirts and the ratio of triples over total number of stars is 
enhanced in the post-collapse phase with respect to isolated 
runs starting from a Plummer model. 

The TV dependence of the number of triples at fixed ini- 
tial binary ratio is more difficult to characterize than the 
dependence on the binary ratio at fixed TV. In fact, while in 
the latter case the structural properties of the cluster are 
fixed, by varying TV the core to ha lf mass ratio is chang- 
ing (see IVesperini fc Chernofil Il99l and HTH, THH), so 
that the analysis in terms of Eqs. (J2][3j) is complicated by 
the necessity to integrate over the density profile and by 
the N-dependence of the a coefficients that describe the effi- 
ciency of the different formation/destruction channels. Em- 
pirically we can also note the intrinsic large scatter in the 
measured triples abundances (see the dots representing in- 
dividual runs in Figs. I5I8[1 . so that it is hard to draw a 
firm conclusion. Tentatively, the TV dependence is logarith- 
mic at most, especially at lower N (in fact the core radius 
decreases approximately as l/log(0.1 ■ TV)). However for the 
last three points (N=4096, 8192, 16384) the data are consis- 
tent with a constant value. In fact the decrease in the core 
radius is compensated by a higher central density and by 
an increased number of binaries over singles in the core (see 



Figure 5. Fractional abundance of triples averaged after the core 
contraction for isolated runs starting from a Plummer model with 
10% primordial binaries and different numbers of particles. Sym- 
bols are as in Fig. [4] 
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Figure 6. Like Fig.[5]but for simulations with a tidal field starting 
from a King Wo = 7 model. 



Fig. 21 in HTH), so that the fractional abundance of triples 
should stay approximately constant. 

The evolution proceeds in a qualitatively similar way 
in our runs initialized with an extended binary binding en- 
ergy range (up to 10 4 kT ). The only difference is that, all 
other conditions being equal, triples formation is enhanced 
by a factor between 1.5 and 2 when harder primordial bi- 
naries are present. In fact, binaries that have binding en- 
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Figure 7. Fractional abundance of triples averaged before the 
core contraction for isolated runs starting from a Plummer model 
with 10% primordial binary and different number of particles. 
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Table 1. Asymptotic triple fraction in our equal mass runs. The 
first column reports the number of particles N, the second the 
initial density profile (Isolated Plummer model - PI or tidally 
limited king model - Wo), the binary fraction / is in the third 
column. The fourth and fifth entries are the average asymptotic 
triple fraction and its standard deviation, computed over N rung 
realizations (last entry). 
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Figure 8. Like Fig.[7]but for simulations with a tidal field starting 
from a King Wo = 7 model. 

ergies above as 700 kT are dynamically inert and behave 
essentially like singles during gravitational interactions with 
other stars. Therefore the encounter of one of these ultra- 
hard binaries with a regular (hard) binary can lead to an 
exchange encounter where one component of the regular bi- 
nary is replaced by the ultra-hard pair, much like during a 
single-regular binary interaction. A more efficient produc- 
tion channel for stable triples is therefore available in these 
simulations. Therefore neutron star binaries with a degen- 



erate companion and with short orbital periods have an en- 
hanced probability of ending up being in a triple system 
compared to main sequence binaries. 

To summarize, in a typical globular cluster (N = 3 TO 5 ) 
we expect an order of at least 100 triples if we conservatively 
assume a binary fraction of 10% (with a standard binding 
energy range, up to contact main sequence binaries) and 
Nt/Nj, ~ 5 ■ 10~ 3 (this accounts for the effect of a mass spec- 
trum and tidal field in setting Nr/Ny, see Fig. [3}. The num- 
ber approximately doubles by adopting an extended binding 
energy range that includes binding energies reached by neu- 
tron star binaries with a degenerate companion. 
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Figure 9. Internal (e; nt ) and external (e ex t) eccentricity distri- 
bution for triples measured in a series of simulations with a tidal 
field starting from King Wo = 3,7,11 models with 20% primor- 
dial binaries. 



Figure 10. Internal (dotted) and external (solid) period distri- 
bution - in units of the period for a lkT binary (po) - for triples 
measured in a series of simulations with a tidal field starting 
from King Wo = 3,7,11 models with 20% primordial binaries 
and N = 16384. 



7 PROPERTIES OF NEWLY FORMED 
TRIPLES 

While for the majority of the standard simulations in our 
sample we recorded only the number of stable triples, we 
have a small number of high resolution (N= 16384) simula- 
tions with complete information on the orbital properties 
of triples, saved every 10 td- In addition, all the extended 
binding energy range simulations have recorded information 
of the orbital properties of multiple systems. 

For a sample of 3 simulations with N= 16384 and / = 
20%, Fig. [9] shows the distribution of internal and exter- 
nal eccentricities. The inner eccentricity has a distribution 
peaked at high eccentricities which derives (i) from the in- 
put thermal distributi on of the ecc entricity for primordial 
binaries and (ii) from iKozail l|l962i ) perturbations induced 
by the outer body, which lead to growth of the inner binary 
eccentricity l|Aarsethll200ll ; lAarseth fc Mardlin3l2000l ). Note 
that our simulations are limited to Newtonian dynamics, so 
we are missing relativistic corrections, w hich may alter the 
orbits of the inner binary l|Aarsethl 2007). The outer compo- 
nent of a stable triple is instead biased toward more circular 
orbits. This means that there is a selective disruption of 
triple systems with eccentric outer orbits. 

For the same simulations the period distribution is re- 
produced in Fig. [10] Here we can note that the internal 
component has an orbital period between two and three or- 
ders of magnitude shorter than the outer component. The 
internal orbit corresponds on average to a binary with bind- 
ing energy of a few hundreds of kT while the external orbit 
has a binding energy of a few kT. The distribution of orbital 
periods has a lower limit at log (p/po) = —4.25, where po is 
the orbital period of a lkT binary. This lower limit corre- 



sponds to a binding energy of ~ 700fcT, i.e. a quasi-contact 
binary with IMq main sequence components. 

The dynamical interactions in the star cluster lead to 
a fast destruction of triples with external binding energy 
below lkT; the destruction of the typical triples found in our 
simulations (with external binding energy of 5kT) happens 
on a timescale of about 10 half mass relaxation times. 

The eccentricity distribution in our extended kT range 
runs is similar to that plotted in Fig. [9] showing only a 
marginally less marked circularization of the outer orbit. 
The distribution of the orbital periods is instead different 
(see Fig. Hip . Stable triples are composed preferentially of 
an inner ultra hard binary (Et > 500 kT) and an external 
component that, on average, has a period only marginally 
shorter than in the standard runs. 

The limited statistics in our simulations do not allow us 
to empirically constrain the N-dependence of the dynamical 
properties of triples. On the basis of theoretical modeling we 
do not expect significant variations with N for single mass 
star clusters, so that the orbital properties discussed here 
should be representative for simulations of globular clusters 
with a realistic number of particles. However one impor- 
tant caveat must be made as there is no guarantee that 
these properties are representative for the triple population 
formed in a realistic model that includes a mass spectrum 
and stellar evolution. 



8 PRIMORDIAL TRIPLES 

If a star cluster originates with a (small) population of pri- 
mordial triples, how does the number of surviving primordial 
triples after several billions of years compare with the num- 
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Figure 11. Internal (dotted) and external (solid) period distribu- 
tion - like in Fig, 1101 - for triples measured in a series of simulations 
with a tidal field starting from King Wo = 7 models with 20% 
primordial binaries, N = 8192 and an extended initial binding 
energy range ([5 kT : 10 4 kT\), representative for neutron star 
binaries with a degenerate companion. 



ber of dynamically formed triples? To answer this question 
we have performed a few simulations with N = 4096 start- 
ing from an isolated Plummer model and with a number of 
primordial triples (ftrip = 5% and ftrip = 10%, where ftrip 
is the initial number fraction of triples). These runs are, to 
the best of our knowledge, the first attempt to investigate 
the dynamical evolution of a star cluster with a number of 
particles and of primor dial triples greater than th e few hun- 
dreds particles used by Ivan den Berk et alj ([2006) . 

We started these runs by initializing the inner binary 
components with the same procedure adopted in Heggie, 
Trenti & Hut (2006), i.e. in an energy range [5 : 680] kT. An 
external component has then been added using the triples 
initialization subroutine within NBODY6 (Aarseth 2003), 
with an energy in the range [0.1 : 100] kT when the external 
energy is defined as 2m/a ext with m being the single particle 
mass and a ext being the external semi-axis. This means that 
the softest external component has a semi-axis 100 times 
larger than the softest inner binary. The evolution of the 
number of triples is shown in Fig. 1121 after a slow start 
(due to the initial large core, and therefore low density of 
triples in the core) primordial triples are steadily burned in 
the first 15 t r h(0) with a rate approximately proportional to 
the number of remaining triples (note that if multiplied by 
a factor 2 the curve for ftrip = 5% reproduces closely the 
behavior of the curve for ftrip = 10%). This agrees well with 
the expectation based on Eq. [3] 

At later times the destruction rate slows down slightly 
in Fig.[l2]due to the expansion of the half mass radius, which 
increases the instantaneous relaxation time with respect to 
the initial one used in the figure. At about 25 t r h(0), a time 
larger than the typical age of a galactic globular cluster, the 



Figure 12. Evolution of the fraction of triples ftrip in two runs 
with N = 4096 starting from an isolated Plummer model and 
ftrip = 5%, 10%. The destruction rate is qualitatively well de- 
scribed by Eq. \3\ 

number of primordial triples is higher than the number of 
dynamically formed triples from 10 — 20% primordial bina- 
ries runs, if the initial fraction of triples is higher than about 
one percent. If we add also primordial binaries in runs with 
primordial triples, the destruction rate of triples is slightly 
enhanced, as expected from Eq.|3](we have verified this with 
a N=1024 run with ftrip = 3% and 7% primordial binaries). 



9 DETECTION OF A TRIPLE SYSTEM WITH 
ONE PULSAR 

As described in Sec. 4, for typical cluster parameters and 
IMq stars, a lkT binary has a semi-major axis of ~10 AU, 
and an orbital period of po~20yrs. Such an orbit is on the 
high-end of the external orbit distribution as shown in Fig- 
ure 9. For an MSP in such a lkT orbit, the instantaneous 
orbitally induced period derivative would be Porbit^W 14 
(compared to a typical intrinsic Ppsr~10 -20 , and 1-yr tim- 
ing accuracies of AP~10 -22 ), and would be identified within 
a couple months of timing observation^. Shorter orbital pe- 
riods would cause significantly larger period derivatives since 
Porbu <x P~ 4 ^ 3 , where p is the orbital period. There is no ev- 
idence for such orbitally induced accelerations in any of the 
~60 globular cluster binary pulsars that have been or are 
currently being timed. 

2 This time is longer than what would be estimated by simply 
considering the six orders of magnitude difference in the accu- 
racy of the 1-yr timing vs. the required accuracy for detection of 
a triple system. This is because there are many degeneracies (for 
example a positional offset) that can mask as a period deriva- 
tive and that require a few months of observations before being 
excluded. 
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Over the past several years, intensive globular cluster 
pulsar searches have uncovered numerous systems that are 
almost certainly due to exchange interactions. These sys- 
tems include at least 10 recycled pulsars in highly eccen- 
tric ( e > 0.25) orbits with periods betwe en ~1— 30 days 
(e.g. iFreire et ai1l2004l ; iRansom et al.ll2005l ). Such systems 
are very similar to the external orbits found in the triples 
described in Figure 9. However, precise timing observa- 
tions rule out the existence of internal orbits composed of 
at least one main-sequence star to a high degree of con- 
fidence. Surveys have also uncovered another class of ex- 
change products, the 3 pulsar— "main-seq uence" binaries 
NGC6397 A, Terzan 5 P, and Terzan 5 ad dD'Amico et all 
l200ll ; IRansom et ail 120051 ; iHessels et all 120061 ). These svs- 
tems have circular orbits of duration 0.3— 1.5 days and 
strange "bloated" companions that cause irregular eclipses. 
The low-eccentricity orbits were likely produced via tidal cir- 
cularization. The erratic timing and irregular eclipses could 
be caused by an internal orbit, but the small orbital sep- 
arations (2—6 Rq) and mass-function constraints indicate 
likely main-sequence dwarf companions. 

The simulations we present suggest there should be ap- 
proximately 100 times more binaries than triples in a typical 
globular cluster (see Figure 2). Given that we know of ~100 
globular cluster MSP binaries and have already one con- 
firmed triple system, the st range MSP-WD-planet system 
in M4 l|Thorsett et al.|[l999l ). the current data are roughly 
consistent with the simulations. The next generation of radio 
telescopes, like the Square Kilometer Array, should uncover 
hundreds of new globular cluster MSPs, greatly increasing 
our chances of finding a pulsar triple system comprised of 
three stars of comparable mass. 



10 DISCUSSION 

In this paper we have presented a systematic investigation of 
the frequency of dynamically formed triples for an extensive 
set of direct simulations starting with a significant fraction 
of primordial binaries. 

We have shown that on a timescale of a few relaxation 
times the abundance of triples reaches an equilibrium value 
that depends linearly on the primordial density of binaries 
and that is about two orders of magnitude smaller. Sim- 
ulations including a tidal field present a relatively larger 
fraction of triples (N t /N) with respect to isolated runs with 
similar initial conditions. The presence of a mass spectrum 
tends instead to reduce the triple fraction. 

Stable triple stars in our standard runs, representative 
for main sequence stars, are primarily formed due to four 
body encounters that lead to the escape of a single star. 
The formation mechanism influences the properties of these 
systems, which are primarily made of a hard inner binary, 
and which in a typical globular cluster would have an aver- 
age orbital period of a few days, accompanied by an external 
star with a period on average between 100 and 1000 times 
longer. The inner component is also, on average, on a slightly 
more eccentric orbit than the outer member of the triple. In 
fact eccentric outer orbits are preferentially perturbed by 
gravitational encounters or break up the inner binary. 

Neutron stars binaries with a degenerate companion can 
reach binding energy well above the 700 kT limit of two main 



sequence stars. To quantify the rate of formation of stable 
triples involving two degenerate stars we have carried out a 
series of simulations with binding energy distributed up to 
10 4 kT. In these runs we observe an enhanced triple abun- 
dance, as an ultra-hard binary made of two degenerate stars 
essentially behaves like a single body during gravitational 
encounters and can efficiently interact with a standard bi- 
nary to form a triple via an exchange-like encounter. In a 
typical globular cluster, the inner orbital period may in this 
case be of the order of a few hours, while the outer compo- 
nent of the system is expected to be of the order of a few 
years, much like in the case of triples with main sequence 
stars. 

Triple stars with at least one pulsar component are ex- 
pected to have a dynamical origin, as it is highly unlikely 
that a pre-existing triple system can avoid being disrupted 
by the Supernova explosion that forms the pulsar. Therefore 
we expect that the properties of pulsar triples are in agree- 
ment with those derived by runs starting with primordial 
binaries only, where all the triples have been formed trough 
stellar encounters (though a single pulsar could become a 
member of a primordial triple via an exchange process). To 
date, ~100 pulsar binaries are known in the globular cluster 
system (Camilo & Rasio 2005). Based on our idealized sim- 
ulations, we expect that the discovery of a pulsar in a triple 
system of comparable masses is likely to happen soon. How- 
ever, we note that the majority of the known pulsar binaries 
have companions of mass ~10% of the pulsar mass rather 
than ~ 1 Mq. How this difference affects our conclusions is 
currently unknown and will demand more detailed studies 
with a spectrum of stellar masses beyond that used in this 
study. 

Higher order hierarchical systems, such as quadruplets 
and quintuplets are also occasionally observed during our 
runs, and their average number density is about two orders 
of magnitude smaller than that of triples. This makes a 
more detailed characterization of their properties extremely 
challenging, at least for simulations with only a limited 
number of particles like ours. A survey aimed at quantifying 
the observed fraction of multiple stars in globular clusters 
would allow us to understand if the observed frequency 
of triples and higher order systems is consistent with 
dynamical production from primordial binaries. Our inves- 
tigation suggests that a triple to binary ratio up to ~ 2% is 
consistent with a purely dynamical origin. If the observed 
number of triples is higher, then primordial triples have to 
be introduced in the theoretical model ing, as investigated 
for sm all N open clusters (N = 182) by| van den Berk et alJ 
|2006t ). 

Acknowledgements 

We thank the referee for a careful reading of the 
manuscript and John Fregeau for useful suggestions. MT 
was supported in part by NASA award HST-AR.10982. 



REFERENCES 

Aarseth S., 2000, in The Formation of Hierarchical Systems 

in Star Clusters Schielicke, R. E. ed., AGM, 16, T0 6A 
Aarseth, S. & Mardling R. A. 2000, (astrrJph/0011514 



Triple Stars with and without a Pulsar in Star Clusters 11 



Aarseth S., 2001, in Dynamics of Star Clusters and the 
Milky Way, ASP Conference Series, 228, 111 

Aarseth S., 2003, Gravitational N-body Simulations. Cam- 
bridge University Press 

Aarseth S. 2004, in "The environment and evolution of dou- 
ble and multiple stars, ed. Allen & Scarfe, IAU Colloquim 
191. 

Aarseth, S. 2007, MNRAS, 378, 285 

Albrow M. D., Gilliland R. L., Brown T. M., Edmonds 
P. D., Guhathakurta P., Sarajedini A., 2001, Ap.J., 559, 
1060 

Baumgardt, H. and Makino, J. 2003, MNRAS, 340, 227 
Bellazzini M., Fusi Pecci F., Messineo M., Monaco L., Rood 

R. T., 2002, A.J., 123, 1509 
Camilo, F. and Rasio, F.A. 2005, in ASP Conf. Ser. 328: 

Binary Radio Pulsars, ed. F.A. Rasio & I.H. Stairs, 147 
D'Amico, N. and Possenti, A. and Manchester, R. N. and 

Sarkissian, J. and Lyne, A. G. and Camilo, F. 2001, Ap.J., 

561, L89 

de la Fuente Marcos, R. 1996, A&A, 314, 453 

de La Fuente Marcos, R. and Aarseth, S. J. and Kiseleva, 
L. G. and Eggleton, P. P. 1997, Astrophysics and Space 
Science Library, 223, 165 

van den Berk, J., Portegies Zwart, S. F. and McMillan, 
S. L. W. 2006, MNRAS, submitted, |astro-ph/0607456"1 

Freire, P. (2005), in Binary Radio Pulsars, ASP Conference 
Series, F. A. Rasio and I. H. Stairs eds. Vol. 328, 405 

Freire, P. C. and Gupta, Y. and Ransom, S. M. and 
Ishwara-Chandra, C. H. 2004, Ap.J., 606, L53 

Fregeau J. M., Giirkan M. A., Joshi K. J., Rasio F. A., 
2003, Ap.J., 593, 772 

Fregeau, J. M. and Rasio F. A. 2007, Ap.J., 658, 1047 

Fregeau J. M., Giirkan M. A., Rasio F. A., 2005, to ap- 
pear in Few-Body Problem: Theory and Computer Sim- 
ulations, ed C. Flynn, Annales Universitatis Turkuensis; 
also |astro-ph/0512032"1 

Gao B., Goodman J., Cohn H., Murphy B., 1991, Ap.J., 
370, 567 

Giersz M., Spurzem R., 2000, MNRAS, 317, 581 
Giersz M., Spurzem R., 2003, MNRAS, 343, 781 
Heggie D. C, Aarseth S. J., 1992, MNRAS, 257, 513 
Heggie D. C, Mathieu R. D., 1986, in LNP 267: The Use of 
Supercomputers in Stellar Dynamics Standardised Units 
and Time Scales, p 233 
Heggie, D.C & Hut, P. (2003) "The Gravitational Million- 
Body Problem: A Multidisciplinary Approach to Star 
Cluster Dynamics" 2003, Cambridge University Press 
Heggie, D.C, Trenti, M. & Hut, P. (2006), MNRAS, 368, 
677 

Hessels, J. W. T. and Ransom, S. M. and Stairs, I. H. and 
Freire, P. C. C. and Kaspi, V. M. and Camilo, F. 2006, 
Science, 311, 1901 

Hurley, J. R. and Pols, O. R. and Tout, C. A. 2000, MN- 
RAS, 315, 543 

Hurley, J. R. and Pols, O. R. and Aarseth, S. J. and Tout, 

C. A. 2005, MNRAS, 363, 293 
Hut, P. & Trenti, M. 2007, Scholarpedia Encyclopedia of 

Astrophysic^f] 
Kozaim Y. 1962, A.J., 67, 591 



Ivanova, N. and Heinke, C. and Rasio, F. A. and Belczyn- 
ski, K. and Fregeau, J. 2007, MNRAS, submitted, astro- 
ph/0706.4096 

Lorimer, D. R., 2005, Living Rev. Relativity, 8, 7 (Novem- 
ber 9, 2005 version) 

Mardling, R. A. & Aarseth, S. 2001, MNRAS, 321, 3 

McMillan S., Hut P., 1994, Ap.J., 427, 793 

McMillan S., Hut P., Makino J., 1990, Ap.J., 362, 522 

Mikkola, S. 1983, MNRAS, 203, 1107 

Mikkola, S. 1984, MNRAS, 207, 115 

Miller, G. E. and Scalo J. E. 1979, Ap.J.S., 41, 513 

Pulone L., De Marchi G., Covino S., Paresce F., 2003, 
A&A, 399, 121 

Raboud, R. & Mermilliod, J. C. 1998, A&A, 329, 101 

Ransom, S. M. and Hessels, J. W. T. and Stairs, I. H. and 
Freire, P. C. C. and Camilo, F. and Kaspi, V. M. and 
Kaplan, D. L., 2005, Science, 307, 892 

Ransom, S. M., Cordes, J. M., & Eikenberry, S. S. 2003, 
Ap.J., 589, 911 

Spitzer, L. 1987, Dynamical Evolution of Globular Clus- 
ters. Princeton, NJ, Princeton University Press, 1987 

Sterzik, M. F. & Tokovinin, A. A. 2002, A&A, 384, 1030 

Thorsett, S. E. and Arzoumanian, Z. and Camilo, F. and 
Lyne, A. G. 1999, Ap.J., 523, 763 

Trenti, M., Heggie, D.C. & Hut, P. (2007), MNRAS, 374, 
344 

Trenti, M., Ardi, E., Mineshige, S. & Hut, P. (2007), MN- 
RAS 374, 857 
Vesperini, E., & Chernoff, D. F. 1994, Ap.J., 431, 231 

This paper has been typeset from a Tp^X/ LTfjX file prepared 
by the author. 



3 http: //www. scholarpedia. org/article/N-body_Simulat ions 



